clc;
clear;

%% main主脚本
tic;
startTime = cputime;
pde = process_pde_NS_benchmark(1);
h = pde.h
dt = pde.dt
basis_type1 = "P2";
basis_type2 = "P1";

mesh = process_mesh(pde, h);

Nt = ceil((pde.domain(6) - pde.domain(5))/dt);
fs1 = process_fespace(mesh, basis_type1);
fs2 = process_fespace(mesh, basis_type2);

% process_showmesh(mesh, fs1, [1 1 0]);
% process_showmesh(mesh, fs2, [1 1 0]);

A11 = term_U_V(mesh, fs1);
A13 = pde.nu*term_grad_U_grad_V(mesh, fs1);
A14 = term_P_div_V(mesh, fs1, fs2);
O1 = sparse(fs2.Nb, fs2.Nb);
O2 = sparse(fs2.Nb, 1);
A = [1/dt*A11+A13, -A14;
    -A14', O1];
U0 = sparse(2*fs1.Nb, 1);

for m = 0:Nt-1
    A12 = term_Uold_grad_U_V(mesh, fs1, U0);
    AN = blkdiag(A12, O1);
    At = A + AN;
    bt = [1/dt*A11*U0; O2];
    [At, bt] = process_boundary_benchmark(At, bt, mesh, fs1, fs2, pde);
    SOL = At\bt;
    U0 = SOL(1:2*fs1.Nb);
end

process_field(fs1, fs2, SOL);

TOCtime = toc;
CPUtime = cputime - startTime;
disp(["TOC time(s)", num2str(TOCtime)]);
disp(["CPU time(s)", num2str(CPUtime)]);


%% process_pde_NS_benchmark
function pde = process_pde_NS_benchmark(index)
switch index
    case 1
        pde.index = index;
        pde.nu = 1/30;
        pde.u1 = @(x,y,k) (k == 1)*0 + (k == 2)*0 + (k == 3)*1 + (k == 4)*0;
        pde.u2 = @(x,y,k) (k == 1)*0 + (k == 2)*0 + (k == 3)*0 + (k == 4)*0;
        pde.domain = [0, 1, 0, 1, 0, 1];
        pde.bdary = [0, 0, 0, 0];
        pde.h = 1/32;
        pde.dt = 1/100;
    case 2
        pde.index = index;
        pde.nu = 1/2;
        pde.u1 = @(x,y,k) (k == 1)*0 + (k == 2)*0 + (k == 3)*1 + (k == 4)*0;
        pde.u2 = @(x,y,k) (k == 1)*0 + (k == 2)*0 + (k == 3)*0 + (k == 4)*0;
        pde.domain = [0, 1, 0, 1, 0, 1];
        pde.bdary = [0, 0, 0, 0];
        pde.h = 1/40;
        pde.dt = 0.1*pde.h;
    case 3
        pde.index = index;
        pde.nu = 0.0015;
        pde.u1 = @(x,y,k) (k == 1)*0 + (k == 2)*0 + (k == 3)*1 + (k == 4)*0;
        pde.u2 = @(x,y,k) (k == 1)*0 + (k == 2)*0 + (k == 3)*0 + (k == 4)*0;
        pde.domain = [0, 1, 0, 1, 0, 1];
        pde.bdary = [0, 0, 0, 0];
        pde.h = 1/40;
        pde.dt = 0.1*pde.h;
    otherwise
        error("Invalid pde index.");
end
fprintf("pde need h (for example 1/32, 1/64)\n");
fprintf("pde domain: xmin,xmax,ymin,ymax,tmin,tmax\n");
disp(pde.domain);
fprintf("pde boundary of u: (0=Dirichlet,1=Neumann)\n");
disp(pde.bdary);
end

%% process_boundary_benchmark
function [A, b] = process_boundary_benchmark(A, b, mesh, fespace1, fespace2, pde)
Nb1 = fespace1.Nb;
% Dirichlet
for k = [1 2 4 3]
    if 0 == pde.bdary(1,k)
        bdn1 = bdary_node(mesh, fespace1, k);
        A(bdn1,:) = 0;
        A(sub2ind(size(A), bdn1, bdn1)) = 1;
        b(bdn1) = pde.u1(fespace1.node(bdn1,1),fespace1.node(bdn1,2),k);
        bdn2 = bdn1 + Nb1;
        A(bdn2,:) = 0;
        A(sub2ind(size(A), bdn2, bdn2)) = 1;
        b(bdn2) = pde.u2(fespace1.node(bdn1,1),fespace1.node(bdn1,2),k);
    end
end
% stablization
M = term_P_Q(mesh, fespace2);
matrix_A = blkdiag(sparse(2*Nb1,2*Nb1), 1e-6*M);
A = A + matrix_A;
end